/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
#include "HiMusicWaveForm.h"
#include "Rtypes.h"
#include <iostream>
#include <cmath>
#include "TObject.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TF1.h"
#include "TRandom3.h"

TH1F *h;

HiMusicWaveForm::HiMusicWaveForm() : TObject(), fNPoints(4096), fTStep(25), fType(0)
{
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = 0;
}

HiMusicWaveForm::HiMusicWaveForm(Int_t NP, Float_t dt, Float_t data, Bool_t t0) : TObject()
{
    fNPoints = NP;
    fTStep = dt;
    fType = t0;
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = data;
}

HiMusicWaveForm::HiMusicWaveForm(const HiMusicWaveForm &w) : TObject(w)
{
    fNPoints = w.fNPoints;
    fTStep = w.fTStep;
    fType = w.fType;
    // fData=new Float_t[fNPoints];
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = w.fData[i];
}

void HiMusicWaveForm::SetBaseLine(Float_t data)
{
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = data;
}
void HiMusicWaveForm::Reset()
{
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = 0;
}
void HiMusicWaveForm::Print()
{
    for (Int_t i = 0; i < fNPoints; i++)
    {
        std::cout << "Data [" << i << "]:    " << fData[i] << std::endl;
    }
}
void HiMusicWaveForm::Print(Int_t Nstart, Int_t Nstop)
{
    if (Nstart < 0)
        Nstart = 0;
    if (Nstop > fNPoints)
        Nstart = fNPoints;
    for (Int_t i = Nstart; i <= Nstop; i++)
    {
        std::cout << "Data [" << i << "]:    " << fData[i] << std::endl;
    }
}

void HiMusicWaveForm::Plot(Int_t ncolor, TString option, TString type)
{
    if (ncolor == 5)
        ncolor = kYellow + 1;
    TCanvas *can = (TCanvas *)gROOT->FindObject("HiMusicPlotWave");
    if (!can)
        can = new TCanvas("HiMusicPlotWave", "HiMusicPlotWave", 600, 600);
    // can->cd();
    // can->Clear();
    Float_t *time = new Float_t[fNPoints];
    for (Int_t i = 0; i < fNPoints; i++)
    {
        if (type.Contains("point"))
            time[i] = i;
        else
            time[i] = i * fTStep;
    }
    TGraph *gr = new TGraph(fNPoints, time, fData);
    gr->SetTitle("");
    gr->SetLineColor(ncolor);
    gr->SetMarkerColor(ncolor);
    gr->Draw(option);
}
Float_t HiMusicWaveForm::GetData(Int_t i)
{
    if (i < 0 || i >= fNPoints)
    {
        std::cout << "Error! Required Index " << i << " is illegal!" << std::endl;
        return -1;
    }
    else
        return fData[i];
}
void HiMusicWaveForm::SetData(Float_t *data)
{
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = data[i];
}
void HiMusicWaveForm::SetData(std::vector<uint16_t> data)
{
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = data[i];
}
void HiMusicWaveForm::SetData(Int_t i, Float_t data)
{
    if (i < 0 || i >= fNPoints)
    {
        std::cout << "Error! Required Index is illegal!" << std::endl;
        return;
    }
    fData[i] = data;
}
TH1F *HiMusicWaveForm::MakeHist(const char *name, TString flag)
{
    Int_t interval = 1;
    if (flag.EqualTo("time"))
        interval = fTStep;
    h = (TH1F *)gROOT->FindObject(name);
    if (!h)
        h = new TH1F(name, name, fNPoints, 0, fNPoints * interval);
    for (Int_t i = 0; i < fNPoints; i++)
        h->SetBinContent(i + 1, fData[i]);
    return h;
}
TGraph *HiMusicWaveForm::MakeGraph(const char *name, TString type)
{
    Float_t *time = new Float_t[fNPoints];
    for (Int_t i = 0; i < fNPoints; i++)
    {
        if (type.Contains("point"))
            time[i] = i;
        else
            time[i] = i * fTStep;
    }
    TGraph *gr = new TGraph(fNPoints, time, fData);
    gr->SetName(name);
    return gr;
}
TH1F *HiMusicWaveForm::MakeHistWithRange(const char *name, TString flag, Float_t xmin, Float_t xmax)
{
    Int_t interval = 1;
    if (flag.EqualTo("time"))
        interval = fTStep;
    h = (TH1F *)gROOT->FindObject(name);
    if (xmin >= xmax)
    {
        std::cout << "wrong range, xmin and xmax!!!" << std::endl;
        return h;
    }
    Int_t np = (xmax - xmin) / interval;
    // if(!h)
    h = new TH1F(name, name, np, xmin, xmax);
    for (Int_t i = xmin / interval; i < xmax / interval; i++)
    {
        h->SetBinContent(i - xmin / interval + 1, fData[i]);
        // std::cout<<i-xmin/interval+1<<fData[i]<<std::endl;
    }
    return h;
}
void HiMusicWaveForm::Inverse()
{
    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = -fData[i];
}
void HiMusicWaveForm::InverseWithPed(Float_t ped)
{

    for (Int_t i = 0; i < fNPoints; i++)
        fData[i] = 2 * ped - fData[i];
}
void HiMusicWaveForm::Smooth(Int_t n)
{
    if (n % 2 == 0) /// n is a even number
    {
        n = n + 1;
        std::cout << "n is a even number, change it to the next odd number automaticlly." << std::endl;
    }
    if (n % 2 == 1) /// n is a odd number
    {
        Int_t N = n / 2;
        for (Int_t i = 0; i < fNPoints; i++)
        {

            if (i >= N && i < fNPoints - N)
            {
                Float_t tmp = 0;
                for (Int_t j = -N; j <= N; j++)
                {
                    tmp = tmp + fData[i + j];
                }
                tmp = tmp / n;
                SetData(i, tmp);
            }
            else
                SetData(i, fData[i]);
        }
    }
}
Float_t HiMusicWaveForm::GetPedSlope(Int_t Npoints, TString flag)
{
    Int_t interval = 1;
    if (flag.EqualTo("time"))
        interval = fTStep;
    TH1F *hslop = MakeHist("hPedSlop", flag);
    Int_t Nbin = hslop->GetNbinsX();
    Double_t min = hslop->GetBinLowEdge(1);
    Double_t max = hslop->GetBinLowEdge(Nbin);
    TF1 *fpol1 = new TF1("fpol1", "[0]+[1]*x", min, max);
    hslop->Fit(fpol1, "", "", 0, Npoints * interval);
    // hslop->Fit(fpol1,"QN","",0,Npoints*interval);
    Float_t slope = (Float_t)fpol1->GetParameter(1);
    return slope;
}
Float_t HiMusicWaveForm::GetPedestal(Int_t Npoints)
{
    Float_t ped = 0;
    for (Int_t i = 0; i < Npoints; i++)
    {
        ped = ped + fData[i];
    }
    ped = ped / Npoints;

    return ped;
}
void HiMusicWaveForm::Add(HiMusicWaveForm &w)
{
    for (Int_t i = 0; i < fNPoints; i++)
        SetData(i, fData[i] + w.GetData(i));
}
void HiMusicWaveForm::Substract(HiMusicWaveForm &w)
{
    if (fNPoints != w.GetNumOfPoints())
    {
        std::cout << "data points are not equal between two waves.";
        return;
    }
    for (Int_t i = 0; i < fNPoints; i++)
        SetData(i, fData[i] - w.GetData(i));
}
void HiMusicWaveForm::AddPed(Float_t ped)
{
    for (Int_t i = 0; i < fNPoints; i++)
        SetData(i, fData[i] + ped);
}
void HiMusicWaveForm::SubstractPed(Float_t ped)
{
    for (Int_t i = 0; i < fNPoints; i++)
    {
        fData[i] = fData[i] - ped;
    }
}
Float_t HiMusicWaveForm::GetSum(Float_t th)
{
    Float_t sum = 0;
    for (Int_t i = 0; i < fNPoints; i++)
    {
        if (fData[i] > th)
            sum += fData[i];
        // std::cout<<"data["<<i<<"]: "<<fData[i]<<", sum: "<<sum<<std::endl;
    }
    return sum;
}
Float_t HiMusicWaveForm::GetPeak()
{
    Float_t peak = fData[0];
    for (Int_t i = 0; i < fNPoints; i++)
    {
        if (fData[i] > peak)
            peak = fData[i];
    }
    return peak;
}
TH1F *HiMusicWaveForm::GetDiff(Int_t npoint)
{
    Int_t Nbin = fNPoints / npoint;
    // TH1F* hdif=new TH1F("hdif","hdif",Nbin, 0,fNPoints);
    TH1F *hdif = new TH1F("hdif", "hdif", Nbin, 0, fNPoints * fTStep);
    Float_t sum_1 = 0; // sum of previous npoint values
    Float_t sum_2 = 0; // sum of current npoint values
    Float_t diff = 0;
    for (Int_t i = 0; i < Nbin; i++)
    {
        sum_2 = 0;
        for (Int_t j = 0; j < npoint; j++)
            sum_2 += fData[i * npoint + j];
        if (i > 0)
            diff = sum_2 - sum_1;

        // std::cout<<"point: "<<i<<", sum_1: "<<sum_1<<", sum_2: "<<sum_2<<", diff: "<<diff<<std::endl;

        sum_1 = sum_2;
        hdif->SetBinContent(i + 1, diff);
    }
    return hdif;
}
TH1F *HiMusicWaveForm::GetDiff(TH1F *hh)
{
    TH1F *ht = new TH1F("hdif2", "hdif2", fNPoints, 0, fNPoints * fTStep);
    Float_t *data = hh->GetArray();
    Float_t x1 = data[1]; // sum of previous npoint values
    Float_t x2 = 0;       // sum of current npoint values
    Float_t diff = 0;
    for (Int_t i = 2; i <= fNPoints; i++)
    {
        x2 = data[i];
        diff = x2 - x1;
        // std::cout<<"x2: "<<x2<<", x1: "<<x1<<", x2-x1="<<diff<<std::endl;
        ht->SetBinContent(i, diff);
        x1 = x2;
    }
    ht->SetBinContent(1, 0);
    return ht;
}
Float_t HiMusicWaveForm::GetMaximum()
{
    Float_t max = GetData(0);
    for (Int_t i = 1; i < fNPoints; i++)
    {
        if (max < GetData(i))
            max = GetData(i);
    }
    return max;
}
Float_t HiMusicWaveForm::GetMinimum()
{
    Float_t min = GetData(0);
    for (Int_t i = 1; i < fNPoints; i++)
    {
        if (min > GetData(i))
            min = GetData(i);
    }
    return min;
}
Float_t HiMusicWaveForm::GetMinimum(Float_t tlow, Float_t thigh)
{
    Int_t start = std::ceil(tlow / fTStep);
    Int_t stop = std::floor(thigh / fTStep);
    if (start < 0)
        start = 0;
    if (stop >= fNPoints)
        stop = fNPoints - 1;

    Float_t min = GetData(start);
    for (Int_t i = start + 1; i <= stop; i++)
    {
        if (min > GetData(i))
            min = GetData(i);
    }
    return min;
}
TH1F *HiMusicWaveForm::GetSigma(TH1F *hh)
{
    Float_t max = hh->GetMaximum();
    Float_t min = hh->GetMinimum();
    Int_t Nbin = (Int_t)1.2 * (max - min);
    TH1F *ht = new TH1F("hSigma", "hSigma", Nbin, 1.2 * min, 1.2 * max);
    // Float_t *data = hh->GetArray();
    Int_t N = hh->GetNbinsX();
    for (Int_t i = 1; i <= N; i++)
    {
        Float_t x = hh->GetBinContent(i);
        ht->Fill(x);
    }
    return ht;
}
Float_t HiMusicWaveForm::SlowFilter(Int_t fl, Int_t fg, Float_t tau) // 梯形成形
{
    // in 输入波形数据
    // out 输出波形数据
    //  n 数据长度
    // fl
    // fg

    Float_t *out = new Float_t[fNPoints];
    Int_t n = fNPoints;
    // double *v=(double*)in;//中间迭代数组
    // double tau=4680;
    Float_t b1 = exp(-1 / tau); // 衰减相关M=exp(-1/tau)
    Float_t c0 = -(1.0 - b1) * pow(b1, (Float_t)fl) / (1 - pow(b1, (Float_t)fl));
    Float_t c1 = (1.0 - b1);
    Float_t c2 = (1.0 - b1) / (1.0 - pow(b1, (Float_t)fl));
    // std::cout<<"b1 c0 c1 c2:"<<b1<<"\t"<<c0<<"\t"<<c1<<"\t"<<c2<<std::endl;

    // baseline
    Float_t bsum0 = 0;
    for (Int_t y = 0; y < fl; y++)
    {
        bsum0 += fData[y];
    }
    Float_t bsum1 = 0;
    for (Int_t y = fl; y < (fl + fg); y++)
    {
        bsum1 += fData[y];
    }
    Float_t bsum2 = 0;
    for (Int_t y = (fl + fg); y < (2 * fl + fg); y++)
    {
        bsum2 += fData[y];
    }
    Float_t baseline = c0 * bsum0 + c1 * bsum1 + c2 * bsum2;
    // std::cout<<"baseline:"<<baseline<<std::endl;

    Int_t offset = 2 * fl + fg - 1;
    // in[0]->in[2*fl+fg]有2*fl + fg+1个数据点（此值仅说明什么时候可以迭代，不影响算法）
    // 算法仅需2*fl + fg个数据点
    for (Int_t x = offset; x < n; x++)
    {
        Float_t fsum0 = 0;
        for (Int_t y = (x - offset); y < (x - offset + fl); y++)
        {
            fsum0 += fData[y]; // 0->fl求和s0
        }
        Float_t fsumg = 0;
        for (Int_t y = (x - offset + fl); y < (x - offset + fl + fg); y++)
        {
            fsumg += fData[y];
        }
        Float_t fsum1 = 0;
        for (Int_t y = (x - offset + fl + fg); y < (x - offset + 2 * fl + fg); y++)
        {
            fsum1 += fData[y]; // fl+fg->2*fl+fg求和s1
        } // std::cout<<"fsum0:"<<fsum0<<"\t"<<"fsumg:"<<fsumg<<"\t"<<"fsum1:"<<fsum1<<std::endl;
        out[x] = c0 * fsum0 + c1 * fsumg + c2 * fsum1 - baseline; // 平滑、作差、调节增益
    }

    for (Int_t x = 0; x < offset; x++)
    {
        out[x] = out[offset];
    }
    SetData(out);
    // std::cout<<out[180+fl+fg]<<"       ch_tau:"<<tau<<std::endl;
    return out[180 + fl + fg];
}
Float_t HiMusicWaveForm::GaussFilter(Int_t k, Float_t baseline) // S-K滤波器
{
    Float_t *out = new Float_t[fNPoints];
    Int_t n = fNPoints;

    for (Int_t y = 0; y < n; y++)
    {
        fData[y] -= baseline;
    }

    out[0] = 2 * fData[0] / (1 + k + k * k);
    out[1] = ((k + 2 * k * k) * out[0] + 2 * fData[1]) / (1 + k + k * k);
    // out[1]=(2*fData[1] )/(1+k+k*k);

    for (Int_t y = 2; y < n; y++)
    {
        out[y] = ((k + 2 * k * k) * out[y - 1] - k * k * out[y - 2] + 2 * fData[y]) / (1 + k + k * k);
    }

    SetData(out);
    return 1;
}
Float_t HiMusicWaveForm::PzcCircuit(Float_t m_k1, Float_t m_k2, Float_t baseline) // 极零相消
{
    Float_t *out = new Float_t[fNPoints];
    Int_t n = fNPoints;

    for (Int_t y = 0; y < n; y++)
    {
        fData[y] -= baseline;
    }

    out[0] = (1 + m_k1) * fData[0] / (1 + m_k1 + m_k2);

    for (Int_t y = 1; y < n; y++)
    {
        out[y] = ((1 + m_k1) * fData[y] - fData[y - 1] + out[y - 1]) / (1 + m_k1 + m_k2);
    }

    SetData(out);
    return 1;
}
Float_t HiMusicWaveForm::CRCircuit(Float_t m_d, Float_t baseline) // CR高通微分电路
{
    Float_t *out = new Float_t[fNPoints];
    Int_t n = fNPoints;

    for (Int_t y = 0; y < n; y++)
    {
        fData[y] -= baseline;
    }

    out[0] = fData[0];

    for (Int_t y = 1; y < n; y++)
    {
        out[y] = m_d * (fData[y] - fData[y - 1]) + m_d * out[y - 1];
    }

    SetData(out);
    return 1;
}
Float_t HiMusicWaveForm::RCCircuit(Float_t m_d, Float_t baseline) // RC低通积分电路
{
    Float_t *out = new Float_t[fNPoints];
    Int_t n = fNPoints;

    for (Int_t y = 0; y < n; y++)
    {
        fData[y] -= baseline;
    }

    out[0] = fData[0];

    for (Int_t y = 1; y < n; y++)
    {
        out[y] = (1 - m_d) * fData[y] + m_d * out[y - 1];
    }

    SetData(out);
    return 1;
}
Float_t HiMusicWaveForm::GMM_EM(Int_t m_Point, Int_t m_k) // EM算法
{                                                         // 设置数据
    static TRandom3 *rn = new TRandom3();
    rn->SetSeed(0);

    Int_t N = m_Point; // 观测数据数总数N 第j个数据
    Int_t K = m_k;     // 模型数K 第k个模型
    Int_t number = 0;  // 迭代次数
    Float_t y[N];      // 数据
    Float_t Yjk[K][N]; // 数据响应度
    Float_t Uk[K];     // 模型均值响应度
    Float_t Dk[K];     // 模型方差响应度
    Float_t Ck[K];     // 函数权重
    Float_t Q[K];      // 数据函数值

    Float_t ped = 0; // 基线

    // 例：SumCkQ=SUM(Ck[K]*Q[K])
    Float_t fabsDU;        // 中间值 迭代中止条件
    Float_t SumUk = 1;     // 中间值 迭代中止条件
    Float_t SumCkQ = 0;    // 中间值
    Float_t SumYjk[K];     // 中间值
    Float_t SumYjky[K];    // 中间值
    Float_t SumYjky_Uk[K]; // 中间值

    // 读数据
    Float_t u = 0, d = 0;
    for (Int_t j = 0; j < N; j++)
    {
        y[j] = GetData(j);
        u += y[j] / N;
    }
    for (Int_t j = 0; j < N; j++)
    {
        d += ((y[j] - u) * (y[j] - u));
    }
    d = d / N;

    // 设置迭代初值
    for (Int_t k = 0; k < K; k++)
    {
        Uk[k] = y[k];
        Dk[k] = fabs(d);
        if (Dk[k] == 0)
            std::cout << "in_hi\t" << Dk[k] << std::endl;
        Ck[k] = 1. / N;
    }

    while (SumUk > 1E-6 && number < 500)
    {
        number++;
        // E步
        for (Int_t k = 0; k < K; k++)
        { // 中间值重置
            SumCkQ = 0;
            SumYjk[k] = 0;
            SumYjky[k] = 0;
            SumYjky_Uk[k] = 0;
        }
        for (Int_t k = 0; k < K; k++)
        { // 遍历模型K
            // if(Dk[k]==0)std::cout<<"hi,Dk[k]->Q\t"<<Dk[k]<<std::endl;
            for (Int_t j = 0; j < N; j++)
            { // 遍历数据N
                SumCkQ = 0;
                for (Int_t q = 0; q < K; q++)
                { // 遍历模型k求数据对应函数值
                    Q[q] = (1 / sqrt(6.28 * Dk[q])) * exp(-(y[j] - Uk[q]) * (y[j] - Uk[q]) / (2 * Dk[q]));
                    // if(Q[q]==0)std::cout<<"U["<<q<<"]:"<<Uk[q]<<"\tD["<<q<<"]:"<<sqrt(Dk[q])<<"\tC[0]、C[1]"<<Ck[0]<<"  "<<Ck[1]<<"\ty"<<y[j]<<std::endl;
                    SumCkQ += Ck[q] * Q[q]; // if(SumCkQ==0)std::cout<<"hi\t"<<SumCkQ<<std::endl;
                } // if(SumCkQ==0)std::cout<<"hi,SumCkQ==0->Yjk[k][j]:error\t"<<std::endl;
                Yjk[k][j] = Ck[k] * Q[k] / SumCkQ;                            // if(Yjk[k][j]==0)std::cout<<"hi,Yjk[k][j]==0 Maybe,SumYjk[k]==0\t"<<Yjk[k][j]<<std::endl;
                SumYjk[k] += Yjk[k][j];                                       // 中间值
                SumYjky[k] += Yjk[k][j] * y[j];                               // 中间值
                SumYjky_Uk[k] += Yjk[k][j] * (y[j] - Uk[k]) * (y[j] - Uk[k]); // 中间值
            }
            if (SumYjk[k] == 0)
                std::cout << "hi,SumYjk[k]==0->Uk[k]:error\t" << std::endl;
            if (SumYjky_Uk[k] == 0)
                std::cout << "hi,SumYjky_Uk[k]==0->Dk[k]:error\t" << std::endl;
        }

        // M步
        SumUk = 0;
        for (Int_t k = 0; k < K; k++)
        { // 遍历模型K
            fabsDU = Uk[k];

            Uk[k] = SumYjky[k] / SumYjk[k]; // M步

            fabsDU = fabs(fabsDU - Uk[k]);
            SumUk += fabsDU; // 计算迭代中止条件

            Dk[k] = SumYjky_Uk[k] / SumYjk[k]; // M步
            Ck[k] = SumYjk[k] / N;             // M步
        }
        // for(Int_t k=0;k<K;k++){//遍历模型K
        //    if(Ck[k]<1E-6||Dk[k]<1E-3)goto end;
        // }

    } // while_end

    // end://模型退化
    for (Int_t k = 0; k < K; k++)
    { // 遍历模型K 输出结果
        std::cout << "U[" << k << "]:" << Uk[k] << "\tD[" << k << "]:" << sqrt(Dk[k]) << "\tC[" << k << "]:" << Ck[k] << std::endl;
        ped += Uk[k] * Ck[k];
    }
    // std::cout<<"While_number:  "<<number<<"   ped:"<<ped<<std::endl;
    return ped;
}
HiMusicWaveForm HiMusicWaveForm::ReSample(Float_t dt, Int_t Nfrom, Int_t Nto)
{
    if (Nto == -1)
        Nto = fNPoints;
    if ((Int_t(dt) % Int_t(fTStep)) != 0)
        std::cout << " dt should be integer multiples of the fTStep! " << std::endl;
    Int_t Nskip = dt / fTStep; /// points should be skipped during resampling
    HiMusicWaveForm wave;
    wave.SetTStep(dt);
    Int_t newPoints = Int_t((Nto - Nfrom) / Nskip);
    wave.SetNPoints(newPoints);
    // std::cout<<"dt0= "<<fTStep<<", dt= " <<dt<<", newpoints: "<<newPoints<<", Nskip: "<<Nskip<<", sample from: "<<Nfrom <<" to "<<Nto<<std::endl;

    for (Int_t i = Nfrom; i < Nfrom + newPoints; i++)
    {
        wave.SetData(i - Nfrom, fData[(i - Nfrom) * Nskip + Nfrom]);
        // std::cout<<"i: "<<i<<", data0: "<<(i-Nfrom)*Nskip+Nfrom<<" = "<<fData[(i-Nfrom)*Nskip+Nfrom]<<std::endl;
    }
    return wave;
}

ClassImp(HiMusicWaveForm);